!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! the subroutine is used to calcaulate inverse wall 
! papramenters by MOC.
! ref. to chapter 16 of 'Gas Dynamics' by Maurice J. Zucrow & Joe D. Hoffmann
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! coded by : B. G.
! created  : 2015-06-18
! revised  :
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! include  :
!   subroutine : InverseW
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! nomenclature:
!   k13       slope of line 1-3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! warning:
!   1. angle in rad, not degree
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine InverseW()
   use VariableDef
   implicit none
!              4()
!     3
!         (2)
!               1
!     given variable
!   real*8::x1,y1,Vu1,Vv1,x3,y3,Vu3,Vv3,x4,y4,theta4
   integer::icor,i,j
!     temp variable
   real*8::Vu2p,Vv2p,theta2p,V2p,M2p,l2p
   real*8::Vup,Vvp,yp,thetap,Vp,ap,Mp,lp,Qp,Rp,Sp,Tp
   real*8::x2f,y2f,Vu4f,Vv4f,a4
   real*8::k13
!     calculate variable
!   real*8::Vu4,Vv4,V4,p4,rou4,T4,M4
  
   x2f=0.0
   y2f=0.0
   Vu4f=0.0
   Vv4f=0.0
   icor=6
   y1=abs(y1)
   Vv1=abs(Vv1)
   y3=abs(y3)
   Vv3=abs(Vv3)
   y4=abs(y4)
   thetab=abs(thetab)

   Vu2=Vu3
   Vv2=Vv3
   Vu4=Vu3
   Vv4=Vv3
   k13=(y3-y1)/(x3-x1)
   do i=1,icor
!     parameters in point 2
      do j=1,icor
         Vu2p=(Vu2+Vu4)/2
         Vv2p=(Vv2+Vv4)/2
         theta2p=atan(Vv2p/Vu2p)
         V2p=sqrt(Vu2p**2+Vv2p**2)
         call Thermo(V2p,p4,rou4,T4,a4,M2p)
         l2p=tan(theta2p+asin(1/M2p))
         
         x2=(k13*x1-y1+y4-l2p*x4)/(k13-l2p)
         y2=y4+l2p*(x2-x4)
         Vu2=Vu1+(x2-x1)*(Vu3-Vu1)/(x3-x1)
         Vv2=Vv1+(x2-x1)*(Vv3-Vv1)/(x3-x1)
      
!         if (i==1) then
!            Vu4=Vu2
!            Vv4=Vv2
!         end if
         
         if ((abs(x2-x2f) <= 0.000001) .or. (abs(y2-y2f) <= 0.000001)) then
            exit
         end if
         
         x2f=x2
         y2f=y2
      end do
      
!   predictor corfficients
      Vup=(Vu2+Vu4)/2
      Vvp=(Vv2+Vv4)/2
      yp=(y2+y4)/2
      
      thetap=atan(Vvp/Vup)
      Vp=sqrt(Vup**2+Vvp**2)
      call Thermo(Vp,p4,rou4,T4,ap,Mp)
      lp=tan(thetap+asin(1/Mp))
      
      Qp=Vup**2-ap**2
      Rp=2*Vup*Vvp-Qp*lp
      Sp=delta*Vvp*ap**2/yp
      Tp=Sp*(x4-x2)+Qp*Vup+Rp*Vvp
      
      Vu4=Tp/(Qp+Rp*tan(thetab))   ! velocity calcaulated
      Vv4=Vu4*tan(thetab)

      if ((abs(Vu4-Vu4f) <= ev) .and. (abs(Vv4-Vv4f) <= ev)) then
         exit
      end if
      
!      corrector   
      Vu4f=Vu4
      Vv4f=Vv4
   end do
   V4=sqrt(Vu4**2+Vv4**2)
   call Thermo(V4,p4,rou4,T4,a4,M4)

!   for downside wall   
   if (wallflag==-1) then
      y1=-y1
      Vv1=-Vv1
      y3=-y3
      Vv3=-Vv3
      y4=-y4
      thetab=-thetab
      Vv4=-Vv4
   end if

end subroutine InverseW
